////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  ----------------------------- Ebrahim Foulaadvand, 08 July 2012 ------------------------------------- //
//                                                                                                        //
// The programme "DampedOscillator" evaluates the temporal evolution of a damped harmonic oscillator      //
// for the three case of underdamped, overdamped and critically damped by three algorithms:               //
// Euler-Richardson (ER), Euler-Crommer (EC) and 2nd order Runge-Kutta (RK2). When you choose             //
// one of the algorithms, remeber to inactive the corresponding routines/files for the other ones.        //                                                                                               //
//                                                                                                        //
//                                                                                                        //                                                                                                         //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;

double a (double &x,double &v); // a is the accelration function

double  tau=0.01,T=30,ac,k=9.,m=1.,x0=0.,v0=1.,gamma=8.,E,Eudanal,Ecdanal,Eodanal,kx1,kx2,lx1,lx2,tmid,
xmid,vmid,D,A,B,A1,A2,delta,omega0=sqrt(k/m),omega1,omega2,beta=gamma/(2*m);

// m=oscillator mass,k=spring constant, gamma=dmping coefiicien, T= total simulation time, x0=initial position
// v0=initial velocity, omega0=natural frequency, tau=timestep, E=total energy.

int N=(T/tau); // N=number of timesteps

main(){

vector<double> v(N+1,0),x(N+1,0),xanal(N+1,0),vanal(N+1,0),Eud(N+1,0),Ecd(N+1,0),Eod(N+1,0);

// Arrays x and v store the oscillator's position and velocity time step t. Arrays Eud,Ecd and Eod
// store the oscillator total energy at each time step for the three cases of underdamped,
// critically damped and overdamped. Arrays xanal and vanal show the anlytic solutions.

ofstream filexER("xER.plt");    // output file for the position (ER algorithm)
ofstream filevER("vER.plt");    // output file for the velocity (ER algorithm)
ofstream filexEC("xEC.plt");    // output file for the position (EC algorithm)
ofstream filevEC("vEC.plt");    // output file for the velocity (EC algorithm)
ofstream filexRK2("xRK2.plt");  // output file for the position (RK2 algorithm)
ofstream filevRK2("vRK2.plt");  // output file for the velocity (RK2 algorithm)

//ofstream fileEudEC("EudEC gamma=5.plt");
//ofstream fileEcdEC("EcdEC gamma=6.plt");
ofstream fileEodEC("EodEC gamma=8.plt");
//ofstream fileEudanal("Eudanal gamma=5.plt");
//ofstream fileEcdanal("Ecdanal gamma=6.plt");
ofstream fileEodanal("Eodanal gamma=8.plt");


ofstream fileEnergyEC("EnergyEC.plt");
ofstream filephaseEC("PhaseSpaceEC.plt");
ofstream filexudanal("xudanal.plt");
ofstream filexcdanal("xcdanal.plt");
ofstream filexodanal("xodanal.plt");


// -------------------- initial conditions --------------------

x[0]=x0; v[0]=v0;

// -------------------- under damped case ----------------------

//omega1=sqrt(omega0*omega0-beta*beta);

//delta=-atan((v0/x0) + beta)/omega1;

//D=x0/cos(delta);

// -------------------- critically damped case ----------------------

//A=x0; B=v0+beta*x0;

// -------------------- over damped case ----------------------

omega2=sqrt(beta*beta-omega0*omega0);

A1=(x0*(omega2+beta)+v0)/(2*omega2);

A2=(x0*(omega2-beta)-v0)/(2*omega2);



//------------------------ Euler-Cromer Method ------------------------

for(int t=0;t<N;t++){

v[t+1]=v[t] + tau*a(x[t],v[t]);
x[t+1]=x[t] + tau*v[t+1];

E=0.5*m*v[t]*v[t] + 0.5*k*x[t]*x[t];

//Eudanal=0.5*D*D*exp(-gamma*t*tau/m)*( (k+m*beta*beta)*cos(omega1*t*tau+delta)*cos(omega1*t*tau+delta) +
//m*omega1*omega1*sin(omega1*t*tau+delta)*sin(omega1*t*tau+delta) +gamma*omega1*
//sin(omega1*t*tau+delta)*cos(omega1*t*tau+delta) );

//Ecdanal=0.5*exp(-gamma*t*tau/m)*( (k+m*beta*beta)*(A+B*t*tau)*(A+B*t*tau) + m*B*B -2*beta*B*m*(A+B*t*tau) );

Eodanal=0.5*exp(-gamma*t*tau/m)*(A1*A1*exp(2*omega2*t*tau)*(k+m*(omega2*omega2+beta*beta-2*beta*omega2)) +
A2*A2*exp(-2*omega2*t*tau)*(k+m*(omega2*omega2+beta*beta+2*beta*omega2)) +
2*A1*A2*(k-m*(omega2*omega2-beta*beta)) );


filexEC<<t*tau<<" "<<x[t]<<endl;
filevEC<<t*tau<<" "<<v[t]<<endl;
fileEnergyEC<<t*tau<<"  "<<E<<endl;
filephaseEC<<x[t]<<"  "<<v[t]<<endl;

//fileEudEC<<t*tau<<" "<<E<<endl;
//leEcdEC<<t*tau<<" "<<E<<endl;
fileEodEC<<t*tau<<" "<<E<<endl;

//filexudanal<<t*tau<<" "<<D*exp(-beta*t*tau)*cos(omega*t*tau+delta)<<endl;
//filexcdanal<<t*tau<<" "<<(A+B*t*tau)*exp(-beta*t*tau)<<endl;
//filexodanal<<t*tau<<" "<<A1*exp(omega2*t*tau) +A2*exp(-omega2*t*tau)<<endl;

//fileEudanal<<t*tau<<" "<<Eudanal<<endl;
//fileEcdanal<<t*tau<<" "<<Ecdanal<<endl;
fileEodanal<<t*tau<<" "<<Eodanal<<endl;


}//end t


//------------------------ Euler-Richardson Method------------------------

for(int t=0;t<N;t++){

tmid=t*tau + 0.5*tau;
xmid=x[t] + 0.5*v[t]*tau;
vmid=v[t] + 0.5*tau*a(x[t],v[t]);

v[t+1]=v[t] + tau*a(xmid,vmid);
x[t+1]=x[t] + tau*vmid;


filexER<<t*tau<<" "<<x[t]<<endl;
filevER<<t*tau<<" "<<v[t]<<endl;

}//end t


//------------------------- Runge-Kutta 2nd order Method ------------------------

for(int t=0;t<N;t++){

kx1=tau*v[t];
lx1=tau*a(x[t],v[t]);

kx2=tau*(v[t]+lx1);
lx2=tau*a(x[t]+kx1,v[t]+lx1);

x[t+1]=x[t]+0.5*(kx1+kx2);
v[t+1]=v[t]+0.5*(lx1+lx2);

filexRK2<<t*tau<<" "<<x[t]<<endl;
filevRK2<<t*tau<<" "<<v[t]<<endl;

}//end t



cout<<"end."<<endl;getch();

return 0;

}

//----------------------------------------------------------------

double a (double &x,double &v) {


return ac=-k*x/m-(gamma/m)*v;

}